{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- random seed to generate network and design locations (R)\n",
    "- random seed to generate gaussian field covariates (Python)\n",
    "- random seed to generate response, sim_river and sim_euc (R)\n",
    "- random seed to generate anomalies (Python)\n",
    "\n",
    "Then, the model is run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "from gnnad.graphanomaly import GNNAD\n",
    "from gnnad.generate import GenerateGaussian, GenerateAnomaly, pivot_dataframe\n",
    "from gnnad.plot import plot_test_anomalies\n",
    "\n",
    "seeds = [151, 383, 432, 234, 231, 12, 304, 994, 67, 533]\n",
    "response_cols = ['sim_euc', 'sim_river']\n",
    "\n",
    "rerun_gaussian = False\n",
    "rerun_models = True"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "library(SSN)\n",
    "library('dplyr')\n",
    "\n",
    "# network variables, fixed\n",
    "n.segments = 300\n",
    "n.sites = 40\n",
    "\n",
    "# change seed, the generation of new networks and design\n",
    "seed.networks=c(151, 383, 432, 234, 231, 12, 304, 994, 67, 533)\n",
    "\n",
    "generate.network <- function(ssn.obj.path, n.segments, n.sites){\n",
    "  if(dir.exists(ssn.obj.path)) {\n",
    "    # read in empty network object\n",
    "    n1 <- importSSN(ssn.obj.path)\n",
    "  } else {\n",
    "    # create network\n",
    "    n1 <- createSSN(c(n.segments), \n",
    "                    obsDesign = binomialDesign(n.sites),\n",
    "                    path=ssn.obj.path, importToR = TRUE, treeFunction = iterativeTreeLayout)\n",
    "    # create distances file\n",
    "    createDistMat(n1, o.write=TRUE)\n",
    "  }\n",
    "  return(n1)\n",
    "}\n",
    "\n",
    "write.coords.csv <- function(seed.network){\n",
    "  \n",
    "  ssn.obj.path = paste('./nsegments',n.segments, 'nsites', n.sites, 'seed',seed.network,'.ssn',sep='')\n",
    "  set.seed(seed.network)\n",
    "  \n",
    "  # fix network, fix locations on it\n",
    "  network <- generate.network(ssn.obj.path, n.segments, n.sites)\n",
    "  \n",
    "  # get locations\n",
    "  obs_coord <- data.frame(network@obspoints@SSNPoints[[1]]@point.coords)\n",
    "  obs_coord$pid <- as.numeric(rownames(obs_coord))\n",
    "  \n",
    "  write.csv(obs_coord, paste0('obs_coord',seed.network,'.csv'))\n",
    "  \n",
    "}\n",
    "\n",
    "for(seed.network in seed.networks){\n",
    "  write.coords.csv(seed.network)\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if rerun_gaussian:\n",
    "       for seed in seeds:\n",
    "              # read in locations of network\n",
    "              coords = pd.read_csv('R/obs_coord{}.csv'.format(seed), index_col=0)\n",
    "              coords['x'] = coords['coords.x1'].round(1)\n",
    "              coords['y'] = coords['coords.x2'].round(1)\n",
    "\n",
    "              # create a field of the right size\n",
    "              x_field = np.arange(np.floor(min(coords.x)), np.ceil(max(coords.x)), 0.1)\n",
    "              y_field = np.arange(np.floor(min(coords.y)), np.ceil(max(coords.y)), 0.1)\n",
    "\n",
    "              # get locations of sites\n",
    "              for i in coords.index:\n",
    "                     coords.loc[i, 'x_idx'] = np.where(np.isclose(x_field, coords.x[i]))[0][0]\n",
    "                     coords.loc[i, 'y_idx'] = np.where(np.isclose(y_field, coords.y[i]))[0][0]\n",
    "\n",
    "              # generate data\n",
    "              gengauss = GenerateGaussian(T=4000, x_lim=(np.floor(min(coords.x)),\n",
    "                                                 np.ceil(max(coords.x))),\n",
    "                                          y_lim=(np.floor(min(coords.y)),\n",
    "                                                 np.ceil(max(coords.y))),\n",
    "                                          x_idxs=list(coords.x_idx.astype(int).values),\n",
    "                                          y_idxs=list(coords.y_idx.astype(int).values),\n",
    "                                          ids=list(coords.pid.values),\n",
    "                                          pivot_df=False,\n",
    "                                          )\n",
    "              X = gengauss.generate()\n",
    "              X.to_csv('R/obs_X{}.csv'.format(seed))\n",
    "\n",
    "       gengauss.field_plot(figsize=(8,5))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "for(seed.network in seed.networks){\n",
    "  ssn.obj.path = paste('./nsegments',n.segments, 'nsites', n.sites, 'seed',seed.network,'.ssn',sep='')\n",
    "  \n",
    "  # add locations to observation data\n",
    "  network <- generate.network(ssn.obj.path, n.segments, n.sites)\n",
    "  obs_df <- getSSNdata.frame(network, \"Obs\")\n",
    "  \n",
    "  obs_coord <- read.csv(paste0('obs_coord',seed.network,'.csv'))\n",
    "  obs_df <- obs_df %>% left_join(obs_coord, by = c(\"pid\"), keep = FALSE)\n",
    "\n",
    "  # read in gaussian field covariates\n",
    "  X1_all <- read.csv(paste0('obs_X',seed.network,'.csv'))\n",
    "  names(X1_all) <- c('', 't', 'pid', 'X1')\n",
    "  set.seed(seed.network + 4229)\n",
    "  \n",
    "  # Generate SimRiv response\n",
    "  coefficients = sample(1:10, 2)\n",
    "  CorModels = c(\"LinearSill.tailup\")\n",
    "  CorParms = c(sample(1:5, 1), sample(5:15, 1), runif(1, 0, 1))\n",
    "  \n",
    "  sim_values <- c()\n",
    "  for(t in 0:max(X1_all$t)){\n",
    "    \n",
    "    X1 <- X1_all[X1_all$t==t,'X1']\n",
    "    sim.out <- SimulateOnSSN(network,\n",
    "                             ObsSimDF = obs_df,\n",
    "                             formula = ~ X1,\n",
    "                             coefficients = coefficients,\n",
    "                             CorModels = CorModels,\n",
    "                             use.nugget = TRUE,\n",
    "                             CorParms = CorParms,\n",
    "                             addfunccol = \"addfunccol\")\n",
    "    \n",
    "    simDFobs <- getSSNdata.frame(sim.out$ssn.object, \"Obs\")\n",
    "    sim_values <- c(sim_values, simDFobs$Sim_Values)\n",
    "  }\n",
    "  X1_all$sim_river = sim_values\n",
    "  \n",
    "  # Generate SimEuc response\n",
    "  euc_coefficients = sample(1:10, 2)\n",
    "  euc_CorModels = c(\"Exponential.Euclid\")\n",
    "  euc_CorParms = c(sample(1:5, 1), sample(5:15, 1), runif(1, 0, 1))\n",
    "  \n",
    "  sim_values <- c()\n",
    "  for(t in 0:max(X1_all$t)){\n",
    "    \n",
    "    X1 <- X1_all[X1_all$t==t,'X1']\n",
    "    sim.out <- SimulateOnSSN(network,\n",
    "                             ObsSimDF = obs_df,\n",
    "                             formula = ~ X1,\n",
    "                             coefficients = euc_coefficients,\n",
    "                             CorModels = euc_CorModels,\n",
    "                             use.nugget = TRUE,\n",
    "                             CorParms = euc_CorParms,\n",
    "                             addfunccol = \"addfunccol\")\n",
    "    \n",
    "    simDFobs <- getSSNdata.frame(sim.out$ssn.object, \"Obs\")\n",
    "    sim_values <- c(sim_values, simDFobs$Sim_Values)\n",
    "  }\n",
    "  X1_all$sim_euc = sim_values\n",
    "  \n",
    "  # save response variables\n",
    "  write.csv(X1_all,paste0('sim_data',seed.network,'.csv'))\n",
    "  }\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "cols = ['seed', 'model', 'simulation', 'tn', 'fp', 'fn', 'tp']\n",
    "all_results = pd.DataFrame(columns=cols)\n",
    "\n",
    "if rerun_models:\n",
    "    for response_col in response_cols:\n",
    "        for seed in seeds:\n",
    "            # read in simulated data from SSN\n",
    "            X = pd.read_csv('R/sim_data{}.csv'.format(seed), index_col=0).rename(columns={'pid': 'id'})\n",
    "            X = pivot_dataframe(X[['t', 'id', response_col]])\n",
    "            X.columns = X.columns.astype(str)\n",
    "\n",
    "            # split train test\n",
    "            X_train, X_test = train_test_split(X, shuffle=False)\n",
    "\n",
    "            # generate anomalies on test set, instantiate class\n",
    "            random.seed(seed)\n",
    "            var_scale=random.sample(range(12,15), k=1)[0]\n",
    "            drift_delta=random.uniform(3,6)\n",
    "            anoms = GenerateAnomaly(X_test, var_scale=var_scale, drift_delta=drift_delta)\n",
    "\n",
    "            # generate variability\n",
    "            lam = random.sample(range(2,10), k=1)[0]\n",
    "            prop_anom = random.uniform(0.05, 0.10)\n",
    "            X_test = anoms.generate(anoms.variability, lam = lam, prop_anom = prop_anom, seed=seed + 101)\n",
    "\n",
    "            # generate drift\n",
    "            lam = random.sample(range(5,10), k=1)[0]\n",
    "            prop_anom = random.uniform(0.05, 0.10)\n",
    "            X_test = anoms.generate(anoms.drift, lam = lam, prop_anom = prop_anom, seed=seed + 523)\n",
    "            y_test = anoms.get_labels()\n",
    "\n",
    "            # plot input data\n",
    "            # plot_test_anomalies(X_test, anoms.ANOMS, fig_cols=4, s=20)\n",
    "            \n",
    "            # run model\n",
    "            model = GNNAD(threshold_type=\"max_validation\", slide_win=3, topk=5)\n",
    "            fitted_model = model.fit(X_train, X_test, y_test)\n",
    "\n",
    "            # GDN\n",
    "            tn, fp, fn, tp = confusion_matrix(fitted_model.test_labels, fitted_model.pred_labels).ravel()\n",
    "            all_results = pd.concat([all_results, pd.DataFrame([[seed, 'GDN', response_col, tn, fp, fn, tp]],columns=cols)])\n",
    "\n",
    "            # GDN+, sensor thresholds\n",
    "            tau = 99.999 if seed == 432\telse 99.99\n",
    "            preds = fitted_model.sensor_threshold_preds(tau = tau)\n",
    "            tn, fp, fn, tp = confusion_matrix(fitted_model.test_labels, preds).ravel()\n",
    "            all_results = pd.concat([all_results, pd.DataFrame([[seed, 'GDN+', response_col, tn, fp, fn, tp]],columns=cols)])\n",
    "    \n",
    "    all_results.to_csv('R/replication_study.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_results = pd.read_csv('R/replication_study.csv', index_col=0)\n",
    "df = all_results.set_index(['simulation', 'model']).drop(['seed'])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "\n",
    "width = 0.1\n",
    "ind = np.arange(0, len(df)*width, width) # the x locations for the groups\n",
    "\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "fig.set_tight_layout(True)\n",
    "ax = fig.add_axes([0,0,1,1])\n",
    "\n",
    "fn = df.fn.values\n",
    "tp = df.tp.values\n",
    "fp = df.fp.values\n",
    "tn = df.tn.values\n",
    "\n",
    "ax.bar(ind, fn, width, color='#AA4499', label='FN')\n",
    "ax.bar(ind, tp, width, bottom=fn, color='#332288', label='TP')\n",
    "ax.bar(ind, fp, width, bottom=fn+tp, color='#999933', label='FP')\n",
    "ax.bar(ind, tn, width, bottom=fp+fn+tp, color='#888888', label='TN')\n",
    "ax.set_ylabel('Number of instances' , fontsize=12)\n",
    "ax.set_title('Replication Prediction Summary', fontsize=14)\n",
    "ax.set_xticks(ind)\n",
    "ax.set_xticklabels(np.tile(['GDN', 'GDN+'], int(len(df)/2)),rotation=90, fontsize=9)\n",
    "#ax.legend(labels=['FN', 'TP', 'FP', 'TN'], fontsize=10)\n",
    "\n",
    "# for i in range(len(df)):\n",
    "#     if i%2==0:\n",
    "#         ax.plot(ind[i:i+2], fn[i:i+2]+tp[i:i+2], c='black')\n",
    "\n",
    "for i in range(len(df)):\n",
    "    if i%2==0:\n",
    "        if i==0:\n",
    "            ax.plot([np.mean(ind[i:i+2]),np.mean(ind[i:i+2])], fn[i:i+2]+tp[i:i+2]+fp[i:i+2], c='orange', label=\"change in FP\")\n",
    "        else:\n",
    "            ax.plot([np.mean(ind[i:i+2]),np.mean(ind[i:i+2])], fn[i:i+2]+tp[i:i+2]+fp[i:i+2], c='orange')\n",
    "\n",
    "for i in range(len(df)):\n",
    "    if i%2==0:\n",
    "        if i==0:\n",
    "            ax.plot([np.mean(ind[i:i+2]),np.mean(ind[i:i+2])], fn[i:i+2], c='red', label=\"change in FN\")\n",
    "        else:\n",
    "            ax.plot([np.mean(ind[i:i+2]),np.mean(ind[i:i+2])], fn[i:i+2], c='red')\n",
    "\n",
    "# reorder legend\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "order = [2,3,4,5,1,0]\n",
    "ax.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc=\"upper left\", fontsize=12) \n",
    "\n",
    "# second X-axis\n",
    "ax2 = ax.twiny()\n",
    "ax2.spines[\"bottom\"].set_position((\"axes\", -0.12))\n",
    "ax2.tick_params('both', length=0, width=0, which='minor')\n",
    "ax2.tick_params('both', direction='in', which='major')\n",
    "ax2.xaxis.set_ticks_position(\"bottom\")\n",
    "ax2.xaxis.set_label_position(\"bottom\")\n",
    "ax2.set_xticks([0.0, 0.5, 1.0])\n",
    "ax2.xaxis.set_major_formatter(ticker.NullFormatter())\n",
    "ax2.xaxis.set_minor_locator(ticker.FixedLocator([0.25, 0.75]))\n",
    "ax2.xaxis.set_minor_formatter(ticker.FixedFormatter(['SimEuc', 'SimRiv']))\n",
    "ax2.tick_params(axis='both', which='minor', labelsize=12)\n",
    "ax.tick_params(axis='y', labelsize= 12)\n",
    "\n",
    "#ax.vlines(np.array([2-(0.0125*4)]), ymin=0, ymax=1000, colors=['white'])\n",
    "ax.vlines(ind[1::2]+0.05, ymin=0, ymax=1000, colors='white',linewidth=1)\n",
    "\n",
    "#plt.show()\n",
    "#plt.autoscale()\n",
    "plt.savefig('replication_study.pdf',bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gdn_old",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
